%% init
clear all
close all
%syms lambda
syms x y q12_re q12_im q  mu1 mu2 mu3 real


lambda =-1+2i;
lambda0 =  lambda;
rangexyfinal = [-6 6 -6 6];
rangexy = [-8 8 -8 8];
w = 1+3i;
w = w/abs(w)*.9;
%syms w
%% 1.5 cast
r = solve(real(lambda^2) + q + real(w)*q12_re -imag(w)*q12_im -x^2+y^2,...
    abs(lambda)^4 +2*abs(lambda)^2*q - (x^2+y^2)^2 + 2*(real(lambda^2) * (q12_re*real(w)- q12_im*imag(w))- imag(lambda^2)*(imag(w)*q12_re + real(w)*q12_im) ) + (1-abs(w)^2)*(q^2 - (q12_re^2 + q12_im^2)) ,...
    q12_re, q12_im);

r = solve(real(lambda^2) + q + real(w*(q12_re+i*q12_im)) -x^2+y^2,...
    abs(lambda)^4 +2*abs(lambda)^2*q - (x^2+y^2)^2 + ...
    real(lambda)^2*real(w)*q12_re - real(lambda)^2*imag(w)*q12_im - 2*real(lambda)*imag(lambda)*real(w)*q12_im - 2*real(lambda)*imag(lambda)*imag(w)*q12_re - imag(lambda)^2*real(w)*q12_re+imag(lambda)^2*imag(w)*q12_im     ...
    + (1-abs(w)^2)*(q^2 - (q12_re^2 + q12_im^2)) ,...
    q12_re, q12_im);

N=200;
XX = [];
YY = [];
% 0.9... [0.3:2:8.3, 12.3:4:20.3, 30.3:15:60.3, 80.3:40:160.3 200.3 250.3  300.3]
for q0=[0.3:2:8.3, 12.3:4:20.3, 30.3:15:60.3, 80.3:40:160.3 200.3 250.3 300.3] 
    q120_re = subs(r.q12_re,{'q'},{q0});
    q120_im = subs(r.q12_im,{'q'},{ q0});

    for index = 2:length(q120_im)
        close all
        scrsz = get(0,'ScreenSize');
        fig = figure('Position',scrsz);
        str = changedisplay(q120_re(index));
        expr1 = '\*\((.+).*?\)\^';
        [tok mat] = regexp(str, expr1, 'tokens', 'match');

        disp('pod odmocninou')
        if length(tok)
            f1 = inline(char(tok{1}),'x','y');
            [x1,y1,fvalues1]=   implot(f1,rangexy,N,'r');
        end    

       hold on
       str = changedisplay(q0^2 -  (q120_re(index)*q120_re(index) + q120_im(index)*q120_im(index))^(1));
       f3 = inline(char(str(7:end)),'x','y')
       [x3,y3,fvalues3] = implot(f3,rangexy,N,'b');
       hold off
       grid on
       [XXX YYY] = ginput(5);        
       XX(end+1) = median(XXX);
       YY(end+1) = median(YYY);
    end
end
XX = [real(lambda0) XX];
YY = [imag(lambda0) YY];
%save w05.mat XX YY
save( [ 'w_' num2str(abs(w)) '_lambda_' num2str(lambda) '.mat'],'XX','YY')
%% druha cast

load w_0.9_lambda_-1+2i.mat
N=500;
scrsz = get(0,'ScreenSize');
fig = figure('Position',[1 450 450 450]);
a = abs(lambda)^2;
b = real(lambda^2);
c = abs(w)^2;
%f1 = inline(['x^2+y^2-' num2str(abs(lambda)^2) ''],'x','y');
%[x1,y1,fvalues]= implot(f1,rangexy,N,'b');
hold on

%f=inline(['(x^2+y^2)^2-2*' num2str(a) '*(x^2-y^2) +' num2str(a^2-4*imag(lambda)^2*a ) ''],'x','y');
%[x,y,fvalues]= implot(f,rangexy,N,'b');



L= x^2 + mu1*(real(lambda^2) + q + real(w)*q12_re -imag(w)*q12_im -x^2+y^2) + mu2* ...
(    abs(lambda)^4 +2*abs(lambda)^2*q - (x^2+y^2)^2 + 2*(real(lambda^2) * (q12_re*real(w)- q12_im*imag(w))- imag(lambda^2)*(imag(w)*q12_re + real(w)*q12_im) ) ) +...
mu3*(q^2-(q12_re^2+q12_im^2));

rr = solve(diff(L,mu1), diff(L,mu2), diff(L,mu3), diff(L,x), diff(L,q), diff(L,q12_re), diff(L,q12_im), x, q, q12_re, q12_im, mu1, mu2, mu3);
%figure(1);
%clf
hold on
%rangexy = [-14 14 -14 14]; N=100;
%for i=1:length(rr.x)
str = changedisplay(rr.x(end-3));
f = inline(['x+' str(7:end)], 'x','y');
[xr,yr,fvalues] = implot(f,rangexy,N,'b');
%end
indexy = find(yr<=imag(lambda0));
xr = xr(indexy);
yr = yr(indexy);
indexy = find(yr>=-imag(lambda0));
xr = xr(indexy);
yr = yr(indexy);
clf
hold on


indexy = find(xr<=real(lambda0));
x0 = xr(indexy);
y0 = yr(indexy);
xy = [x0; y0];
xy=sortrows(xy',2)';


hold on
% fill([ xy(1,:),xy1(1,:)] ,[ xy(2,:) xy1(2,:) ],[0.9 0.9 0.9],'LineWidth',0.001,'LineStyle','none') 
% [x1,y1,fvalues]= implot(f1,rangexy,N,'b');
hold on

plot(xr,yr,'b','LineWidth',2)

fu = inline('x^2*y^2-(x^2-y^2)*30-94','x','y');
[xu , yu, fuv]=implot(fu,rangexy,N,'r');% grid on; hold on
izap = find(xu <= real(lambda));
xu = xu(izap);
yu = yu(izap);
xuk = xu(find(yu>=0 & yu<10));
yuk = yu(find(yu>=0 & yu<10));
xuz = xu(find(yu<=0));
yuz = yu(find(yu<=0));
clf; hold on
fill([xuz, xy(1,:), fliplr(xuk)] ,[ yuz xy(2,:) fliplr(yuk)],[0.75 0.75 0.75],'LineWidth',0.001,'LineStyle','none') % fill -final
[xu , yu, fuv]=implot(fu,rangexy,N,'r');% vykresleni cerveneho omezeni - final

figure(2);
str = changedisplay(rr.x(7));%end-3
fx1 = inline(['x+' str(7:end)], 'x','y');
[xx1,yy1,fvalues] = implot(fx1,rangexy,N,'b');
indexyxx1 = find(yy1 <= imag(lambda) & yy1 >= -imag(lambda));
xx1 = xx1(indexyxx1);
yy1 = yy1(indexyxx1);
hold on
plot(xx1,yy1,'r','LineWidth',2)

str = changedisplay(rr.x(8));%end-3
fx1 = inline(['x+' str(7:end)], 'x','y');
[xx2,yy2,fvalues] = implot(fx1,rangexy,N,'b');
indexyxx1 = find(yy2 <= imag(lambda) & yy2 >= -imag(lambda));
xx2 = xx2(indexyxx1);
yy2 = yy2(indexyxx1);
hold on
plot(xx2,yy2,'r','LineWidth',2);


% horni cast krivek



str = changedisplay(rr.x(9));%end-3
fx1 = inline(['x+' str(7:end)], 'x','y');
[xx1u,yy1u,fvalues] = implot(fx1,rangexy,N,'b');
indexyxx1 = find(yy1u >= imag(lambda) | yy1u <= -imag(lambda));
xx1u = xx1u(indexyxx1);
yy1u = yy1u(indexyxx1);

indexynula = find(abs(xx1u)<1e-10);
xx1u_kladne = xx1u(find( yy1u >=0 & yy1u <=min(abs(yy1u(indexynula)))));
xx1u_zaporne = xx1u(find( yy1u <=0 & yy1u>=-min(abs(yy1u(indexynula)))));
yy1u_kladne = yy1u(find( yy1u >=0 & yy1u <=min(abs(yy1u(indexynula)))));
yy1u_zaporne = yy1u(find( yy1u <=0 & yy1u>=-min(abs(yy1u(indexynula)))));
hold on
plot(xx1u_kladne,yy1u_kladne,'r','LineWidth',2)
plot(xx1u_zaporne,yy1u_zaporne,'r','LineWidth',2)


str = changedisplay(rr.x(10));%end-3
fx1 = inline(['x+' str(7:end)], 'x','y');
[xx2u,yy2u,fvalues] = implot(fx1,rangexy,N,'b');
indexyxx1 = find(yy2u >= imag(lambda) | yy2u <= -imag(lambda));
xx2u = xx2u(indexyxx1);
yy2u = yy2u(indexyxx1);

indexynula = find(abs(xx2u)<1e-10);
xx2u_kladne = xx2u(find( yy2u >=0 & yy2u <=min(abs(yy2u(indexynula)))));
xx2u_zaporne = xx2u(find( yy2u <=0 & yy2u>=-min(abs(yy2u(indexynula)))));
yy2u_kladne = yy2u(find( yy2u >=0 & yy2u <=min(abs(yy2u(indexynula)))));
yy2u_zaporne = yy2u(find( yy2u <=0 & yy2u>=-min(abs(yy2u(indexynula)))));
hold on
plot(xx2u_kladne,yy2u_kladne,'r','LineWidth',2)
plot(xx2u_zaporne,yy2u_zaporne,'r','LineWidth',2)


figure(1);
%vykresleni modrych krivek - final
plot([fliplr(xx1) fliplr(xx1u_kladne) xx2u_kladne xx2 xx2u_zaporne fliplr(xx1u_zaporne) fliplr(xx1)],[fliplr(yy1) fliplr(yy1u_kladne) yy2u_kladne yy2 yy2u_zaporne fliplr(yy1u_zaporne) fliplr(yy1)],'b','LineWidth',2)


arrow([0 0],[real(lambda0) imag(lambda0)])
arrow([0 0],[real(lambda0) -imag(lambda0)])
arrow([0 rangexyfinal(3)],[0 rangexyfinal(4)],'BaseAngle',45)
arrow([ rangexyfinal(1) 0],[ rangexyfinal(2) 0],'BaseAngle',45)

axis(rangexyfinal);

% vykresleni os a popisku

h = text;
set(h, 'interpreter', 'latex')
if real(lambda)==-1
    text([real(lambda)],[ -imag(lambda)-.5],'$\overline{\lambda}$','interpreter','latex','FontSize',16);
    text([real(lambda)],[ imag(lambda)+.5],'$\lambda$','interpreter','latex','FontSize',16)
elseif real(lambda)==-2
    text([real(lambda)-.25],[ -imag(lambda)-.75],'$\overline{\lambda}$','interpreter','latex','FontSize',16);
    text([real(lambda)-.25],[ imag(lambda)+.75],'$\lambda$','interpreter','latex','FontSize',16)
elseif real(lambda)==3
    text([real(lambda)+.5],[ -imag(lambda)],'$\overline{\lambda}$','interpreter','latex','FontSize',16);
    text([real(lambda)+.5],[ imag(lambda)],'$\lambda$','interpreter','latex','FontSize',16)

else
    text([real(lambda)+.1],[ -imag(lambda)-.5],'$\overline{\lambda}$','interpreter','latex','FontSize',16);
    text([real(lambda)+0.1],[ imag(lambda)+.5],'$\lambda$','interpreter','latex','FontSize',16)
end
text(rangexyfinal(2)-1.25,.5,'R e$s$','interpreter','latex','FontSize',16)
text(.25,rangexyfinal(4)-.5,'$\mbox{Im}s$','interpreter','latex','FontSize',16)
box off
axis off


